# include <stdio.h>
# include <malloc.h>
# include <stdlib.h>

extern double interpol_1d(FILE *lptr, double *x, double *u, int n, double xp, double *du/*, float shift*/)
 { //n : Number of points.
  int i, m, ns = 1;
  double den, dif, dift, ho, hp, w, up;
  double *c, *d, *_x, _du;

  if(lptr != (FILE *)0) fprintf(lptr, "Entering into 'interpol_1d()'...\n");

  _x = x;

  x = (double *)malloc(sizeof(double)*n);
/* Copying the value because the memory pointed to "x" will be distroyed during course of execution. */
  for(i = 0; i < n; i++) x[i] = _x[i];// + shift*(_x[i+1] - _x[i]);

  x--;
  u--;

  dif = xp - x[1];
  dif *= (dif < 0) ? -1:1;

  c = (double *)malloc(sizeof(double)*(n+1));
  d = (double *)malloc(sizeof(double)*(n+1));

  for(i = 1; i <= n; i++)
       {
        dift = xp - x[i];
        dift *= (dift < 0) ? -1:1;

        if(dift < dif){
                       ns = i;
                       dif = dift;
                      }

        c[i] = u[i];
        d[i] = u[i];
       }

  up = u[ns--];

  for(m = 1; m < n; m++)
     {
      for(i = 1; i <= n-m; i++)
          {
           ho = x[i] - xp;
           hp = x[i+m] - xp;
           w = c[i+1] - d[i];

           if((den = ho - hp) == (double)0)
              {
               if(lptr != (FILE *)0) fprintf(lptr, "\nError:: In 'interpol_1d()': Two or more input points are very close.\n");
               exit(-1);
              }

           den = w/den;
           d[i] = hp*den;
           c[i] = ho*den;
          }

      up += ( _du = (2*ns < (n-m) ? c[ns+1]:d[ns--]));
     }

   if(du != (double *)0) *du = _du;

   free(x+1); x = _x;
   free(c);
   free(d);

   if(lptr != (FILE *)0) fprintf(lptr, "...'interpol_1d()' ends.\n");

   return up;
  }
